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SUMMARY 

Several explicit integration algorithms with self-adaptive time integration strategies are developed 
and investigated for efficiency and accuracy. These algorithms involve the Runge-Kutta second order, the 
lower Runge-Kutta method of orders one and two, and the exponential integration method. The algorithms are 
applied to viscoplastic models put forth by Freed and Verrilli and Bodner and Partom for thermal/mechanical 
loadings (including tensile, relaxation, and cyclic loadings). 

The large amount of computations performed showed that, for comparable accuracy, the efficiency of 
an integration algorithm depends significantly on the type of application (loading). However, in general, for 
the aforementioned loadings and viscoplastic models, the exponential integration algorithm with the proposed 
self-adaptive time integration strategy worked more (or comparably) efficiently and accurately than the other 
integration algorithms. Using this strategy for integrating viscoplastic models may lead to considerable saving 
in computer time (better efficiency) without adversely affecting the accuracy of the results. This conclusion 
should encourage the utilization of viscoplastic models in the stress analysis and design of structural 
components. 


INTRODUCTION 

The need for an accurate assessment of the time-dependent, inelastic response of structural 
components has led to the development of numerous constitutive material models, called viscoplastic models. 
Accurate assessment of inelastic response is an essential input for predicting the service life of components. 
Viscoplastic models incorporate all aspects of inelastic deformation, including plasticity, creep, and 
relaxation. These models do not endeavor to artificially split the inelastic strain into plastic and creep 
components. Alternatively, they treat all the inelastic strain (including plasticity, creep, and relaxation) as a 
single time-dependent entity. This “unified” representation of time-dependent, inelastic strain enables these 
models to automatically include the observed interactions at high temperatures among plasticity, creep, and 
relaxation. These models, therefore, provide a more realistic and accurate description and response of the 
time-dependent, inelastic behavior of materials. The input of inelastic strain, calculated by using these 
viscoplastic models, into life prediction methods thus leads to a more accurate determination of the service 
life of structural components. 

The mathematical framework of most viscoplastic models consists of a flow law that relates the 
inelastic strain rate to the deviatoric stress and evolution equations for internal state variables that represent 
the resistance of the material to further inelastic flow. Most models use two internal state variables, a tensor 
and a scalar. The tensorial and scalar variables represent the kinematic and isotropic hardenings of the 
material, respectively. The form of the evolution laws for the internal state variables is based on the well- 
known Bailey-Orowan theory (ref. 1), which states that the high-temperature deformation of a material takes 
place under the influence of two competing mechanisms. One mechanism represents the hardening process 
with accumulated deformation, and the other represents a thermal recovery or softening process proceeding 
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with time at temperature. The evolution rate of internal state variables is then defined as the difference 
between the rates of the hardening and recovery processes. Under steady-state conditions these two 
mechanisms balance each other. 

There is a significant incentive to make viscoplastic models more realistic by incorporating into them 
as much material science as possible. Consequently, the mathematical structure of viscoplastic models is very 
complex. The constitutive equations governing the flow behavior and the evolution of internal state variables 
are, in general, highly nonlinear and mathematically “stiff.” This inherent stiff character of the constitutive 
equations causes major problems during their numerical integration. Present-day structural engineering 
problems, such as those encountered in the aerospace and nuclear industries, involve complex thermal/ 
mechanical loading histories. For the solution of these problems the entire loading history has to be traced, 
requiring these equations to be integrated many times. The cost and computer time involved prohibit using the 
conventional method of employing small time steps for accurate integration of these stiff equations. There is, 
therefore, a strong need to develop efficient and accurate integration algorithms with a self-adaptive, time-step 
strategy that can achieve the desired accuracy and stability over the entire integration range (loading 
histories). 

Keeping the aforementioned requirements in mind, numerous researchers have proposed a number of 
integration algorithms with self-adaptive, time-step control. These algorithms are based on explicit or implicit 
integration methods. Kumar et al. (ref. 2) and Banthia and Mukherjee (ref. 3) present several such integration 
algorithms and show that for the viscoplastic model due to Hart (ref. 4) using a simple, one-step, explicit 
Euler integration scheme with automatic time-step control is advantageous. Gear (ref. 5), Cormeau (ref. 6), 
Hughes and Taylor (ref. 7), and Homberger and Stamm (ref. 8) have suggested special implicit integration 
schemes. Numerous other attempts to integrate these stiff equations have also been made, and useful 
information in this regard is available in the literature (refs. 9 and 10). Walker (refs. 11 and 12) and Freed 
et al. (ref. 13) have proposed some explicit and implicit asymptotic exponential integration algorithms and 
applied them to the integration of viscoplastic models. (Note that the references mentioned in this paragraph 
are by no means exhaustive. They are rather indicative of the importance assigned by various researchers to 
the development of efficient and accurate time integration strategies for use with viscoplastic models.) 

We have been involved with the development of suitable integration strategies for time integration of 
viscoplastic models. A previous work by Arya et al. (ref. 14) reports that a simple, self-adaptive time 
integration strategy involving the explicit forward Euler method, although only conditionally stable, works very 
successfully if the time step is properly controlled at each integration step. The main advantage of using an 
explicit integration scheme is the ease with which it can be implemented. Another favorable aspect of using 
explicit integration schemes is that these do not require the assembly and inversion of Jacobian matrices. This 
fact is of great advantage and significance especially when the number of equations involved is large. The 
advantage of using implicit integration schemes is that they may allow large time steps without affecting the 
stability of the integration method. Cormeau (ref. 6) has, however, shown that several implicit methods suffer 
from the same time-step restrictions as does an explicit method, and in such circumstances the implicit 
integration methods offer no special advantages over the simpler explicit methods. 

In light of the aforementioned facts, several explicit time integration strategies are presented in this 
paper. Detailed investigations were made to judge their efficiency, accuracy, and suitability for application to 
viscoplastic models. The integration strategies were developed and investigated in conjunction with the 
following four explicit integration methods: 

(1) Runge-Kutta, second-order method (RK2 method) 

(2) Lower Runge-Kutta method of order one (or Euler-Cauchy Method) (RK1L method) 

(3) Lower Runge-Kutta method of order two (RK2L method) 

(4) Exponential integration method (in conjunction with Runge-Kutta second-order method) 

(RK2EXP method) 
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Method 1 is currently being used by the author for numerically integrating stiff equations of viscoplastic 
models. Methods 2 to 4 have shown promise for integrating stiff equations (refs. 15 and 16). However, their 
adaptability and applicability for efficiently and accurately integrating the constitutive equations of 
viscoplastic models is yet to be explored. The present study endeavored to do this by developing suitable time 
integration strategies that employ these methods. 

The integration strategies presented herein are general in nature and can be applied to any 
viscoplastic model. However, two viscoplastic models are used to illustrate the applicability of these 
integration strategies. One viscoplastic model is put forth by Freed and Verrilli (ref. 17) and the other by 
Bodner and Partom (ref. 18). Although these two viscoplastic models possess identical mathematical 
frameworks, their constitutive equations have several significant differences in their forms and in the functions 
they involve. These differences, detailed later in the paper, assign different degrees of mathematical 
complexity to the models. 

First, the mathematical frameworks of the viscoplastic models put forth by Freed-Verrilli and Bodner- 
Partom are briefly described. The integration strategies developed and employed in this paper are then 
presented in conjunction with the four explicit integration methods listed earlier. Next, procedures to make the 
integration strategies self-adaptive are presented. The integration strategies are then applied to both 
viscoplastic models for three different types of loadings— tensile, relaxation, and cyclic. Finally, numerical 
results are presented and discussed in order to judge the applicability, efficiency, and accuracy of these 
integration strategies in conjunction with the viscoplastic models. 


VISCOPLASTIC MODELS 
Freed-Verrilli Viscoplastic Model 

Following the general mathematical structure adopted by most viscoplastic models and as mentioned 
earlier, the Freed-Verrilli (F-V) viscoplastic model (ref. 17) also incorporates two internal state variables. One 
a tensorial variable called the backstress, accounts for strain-induced or kinematic hardening. The other, called 
the drag stress (strength), represents isotropic hardening and is a scalar variable. The model is based on the 
assumption that there is no explicit coupling between the dynamic and static recovery terms in an evolutionary 
equation for the internal state. This assumption renders a simpler mathematical framework for the model. The 
model also employs a small displacement and a small strain formulation. 

The total strain rate t-- is assumed to be composed of elastic el , inelastic (including plasticity, 

creep, and relaxation) e~, and thermal el, strain rate components and can be written as 

+ e"j + il, where i,j = 1,2,3 (1) 

In equation (1), the symbol e denotes the strain and the dot over the symbol indicates its derivative with 
respect to time t. Superscripts e, v, and t refer to elastic, viscoplastic (inelastic), and thermal components of 
the strain rate. 

For an isotropic material the elastic strain rate il is related to the stress rate 6- by Hooke’s law 


1 + v . v . s 

e y - e ° ij £ a ** 5 y 


( 2 ) 
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where E is Young’s modulus, v is Poisson’s ratio, and 8 z y is the Kronecker delta. The repeated subscripts in 
equation ( 2 ) and elsewhere imply summation over their range. 

The deviatoric stress Sy is defined by 


S ij “ °ij 3 G kkPij 


( 3 ) 


The effective stress 


E.j has the expression 


z,r s u- B » 


(4) 


where By is the backstress. 

Flow law . — Defining the invariant J 2 by 


h=- 2 hjh 


(5) 


the flow law is written as 



( 6 ) 


The function 0, called the thermal diffusivity function, contains the temperature dependence of the model. It is 
defined by 



if 7 > 0.57V 

fit 





if 7 < 0.57V 

m 


(7) 


Q is the activation energy, k is the Boltzmann constant, T is the absolute temperature, and T m is the melting 
point of the material. 

The function Z, called the Zener-Hollomon parameter, is given as 


J AF n , if F< 1 

| A exp[n(F - 1)], if F > 1 


and A and n are material constants. The function F is defined as 
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(9) 


F = 


h 

D 


where D denotes the drag stress. 

Evolution laws . — The differential equations that govern the growth of the internal state variables 
(backstress Bij and drag stress D ) are 


and 


h= H i 


*5 



D = h\ 


ll 

G 


0KG) 


( 10 ) 


( 11 ) 


where 


I 2 = 0Z 


( 12 ) 


In equations (10) and (11), H, h, and L are inelastic material constants. The constant L denotes the limiting 
value of the backstress at kinematic saturation (ref. 17). The recovery function r is defined by using the 
following equations: 


JO. if D = £> 0 
\R(G), if D>D 0 


(13) 


R(G) = 


Jag" -1 , 

I A exp[n(G 


if G<1 
1)]/G, 


if G >1 


(14) 


and 


G = 


L 


S-D 


(15) 


where S and Do are material constants. 

For the theory to be thermodynamically admissible, the dissipativity condition must be satisfied. This 
results in the following constraint on the function r. 
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where 

*2 = ( 17 ) 

The values of elastic and inelastic material constants appearing in the preceding equations are given in 
table I for copper (ref. 17). 

Bodner-Partom Viscoplastic Model 

The viscoplastic model put forth by Bodner (ref. 19) in 1968 is one of the first unified constitutive 
models that does not require an explicit yield function (or loading/unloading conditions) to govern the inelastic 
response of a material. The mechanistic basis for the model development stems from the work of Gilman 
(ref. 20) and Johnston and Gilman (ref. 21). In the Bodner-Partom (B-P) model the inelastic strain rate is 
considered to be a function of the current value of stress a and an internal state variable that represents a 
measure of hardening or resistance to inelastic flow. Although the general mathematical framework of this 
viscoplastic model is akin to that adopted by other viscoplastic models, it has several significant differences 
from them. For example, in this model the inelastic strain rate is expressed as a function of stress a and not as 
a function of overstress, which is defined as the difference of stress a and the backstress. Also, the measure of 
inelastic hardening in this model is taken to be inelastic (plastic) work and not the inelastic strain as is done 
in other viscoplastic models. 

Several modifications to this model have been proposed from time to time. The modifications of the 
model advanced by Bodner and Partom (refs. 18 and 22) included large deformations and isotropic hardening, 
respectively. Stouffer and Bodner (refs. 23 and 24) included the Bauschinger effect associated with load 
reversals and the generally anisotropic nature of inelastic hardening. Bodner (ref. 25) extended the model to 
include evolution of damage. The model used in the present paper (from ref. 18) includes isotropic and 
directional hardening, thermal recovery of hardening, and general temperature dependence of inelastic flow. 
The constitutive equations of the model are given here. 

Flow law . — The total strain rate is again assumed to be the sum of elastic, inelastic (viscoplastic), 
and thermal strain rate components, and the inelastic strain rate is assumed to have the form of the Prandtl- 
Reuss flow law. In symbols 

£ i j=£*+£jj + £-j , where i, j = 1 , 2, 3 (18) 

^ =0 (19) 


where 


S = a — 5 -Om 
tj ij ^ y k* 


( 20 ) 


In these equations the symbols e, S, o, and 5 and the superscripts e, v, and t have the same meaning as 
defined earlier. A repeated index implies summation over its range, and the parameter X is defined later. 
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Kinetic equation . — The inelastic deformations are governed by the kinetic relation 


D p 2 =Dl exp|-[z 2 /(37 2 )j* | 

with 

Z = Z l + Z D ; 


and 




2 11 11 


2 2 


( 21 ) 


( 22 ) 


(23) 


In these equations Z is an internal state variable that represents the resistance of the material to inelastic flow. 
This state variable in the model is assumed to be composed of isotropic hardening Z l and directional hardening 
ZP components. Also, Dq denotes the limiting strain rate in shear and n is a material constant. 

Isotropic-hardening evolution laws. — The growth or evolution of the isotropic-hardening component is 
governed by 


where 


Zj + (XZ 3 -z')w p -A.Z^Z 1 — Z^/zJ 1 

(24) 

a = m 2 (a, - a)w p sin 0 

(25) 

cos-'fv;..^.) or 0 = cos -1 ) 

(26) 

f = Py /^kfikl ’ V ij = 

(27) 

= a ij/^l°kl a ki' u ij~^ij/ ^kl^kl 

(28) 


with 


z\ 0) = Z 0 ; W p = o i} i]j ; W p (0) = 0; <x(0) = 0 (29) 

In these equations Zo is the initial value of the isotropic hardening variable Z 1 ; Z\ denotes the 
limiting (maximum) value of Z ; ; Z 2 is the fully recovered (minimum) value of Z l ; and Z 3 denotes the limiting 
(maximum) value of Z D ■ The parameter 0 is a measure of nonproportionality of loading, and a is a loading- 
history-dependent variable that provides an additional hardening increment in the evolution equation for Z 1 . 
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The coefficient (Xi corresponds to the limiting value of a, and m \ and m 2 denote the hardening rate 
coefficients of and ZP y respectively. The rate of inelastic (plastic) work is denoted by W . The symbols A\ 

and r\ denote the recovery coefficient and the recovery exponent for ZJ. The second-order symmetric tensor py 
represents the directional-hardening state variable. 

Directional-hardening evolution law. — The evolution equation for the directional-hardening variable 
Py has the same general form as that for isotropic hardening. It, however, has a tensorial character. It is written 
as 


Vij - V tj (30) 

with 

Z° =PyK | - / ; Z D ( 0) = 0; p (/ (0) = 0 (31) 

Here m 2 denotes the hardening rate coefficient of Z D . The symbols A 2 and r 2 represent the recovery coefficient 
and the recovery exponent for ZP . In addition to elastic material constants the model has the following 
inelastic material constants: Dq , Zq, Zj, Z 2 , Z 3 , mi,m 2 , aj , A 2 , r 2 . and n. In most cases one can set r\ = 

r 2 , Aj = A 2 , and Zo = Z 2 . The values of all these constants for the material B1900+Hf, taken from reference 
18, are listed in table n. 


NUMERICAL INTEGRATION METHODS 

To develop efficient, accurate, and stable explicit integration strategies for integrating the constitutive 
(differential) equations of the F-V and B-P viscoplastic models, four explicit numerical integration methods 
were considered. The reasons for considering only the explicit integration methods have been explained in the 
introduction. To keep the presentation simple, the mathematical forms of the integration strategies have been 
presented for only one differential equation. Generalization of these integration strategies for a system of 
differential equations is straightforward and presented subsequently. 

Consider the first-order differential equation 

^ = /(*,?) (32) 

dx 


Let the solution, denoted by y, of this differential equation be known at time t. The mathematical formulas for 
the four methods are given here to calculate the solution at time f+/z where h denotes the increment in time t. 


Runge-Kutta, Second-Order Method (or RK2 Method) 

For equation (32) the approximate solution at t+h using the Runge-Kutta method is given by the well- 
known expressions (see, e.g., ref. 26) 


y(t + h) = y + Uk 1 +^) 


(33) 
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with 


k { =hf(t,y)\ f(t + h,y + k x ) 


Lower Order Runge-Kutta Methods (RK2L and RK1L Methods) 

Fehlberg (ref. 15) has proposed several lower order Runge-Kutta formulas for the solution of a system 
of ordinary differential equations and has applied them to heat transfer problems. He has shown for these 
problems that the low-order, Runge-Kutta methods are considerably faster and of equal or better accuracy than 
the conventional explicit integration formulas. Higher order Runge-Kutta formulas offer no special advantages 
because the increase in the integration step size is restricted by the stability considerations resulting from the 
exponential character of the solution. 

It is well established that small truncation errors lead to efficient Runge-Kutta integration formulas 
and that the permissible integration step size is inversely related to the magnitude of these truncation errors. 
The attempt in the lower order integration formulas proposed by Fehlberg has, therefore, been to keep the 
truncation errors as small as possible. In this study we selected two lower order Runge-Kutta methods proposed 
by Fehlberg: the lower Runge-Kutta method of order two (or RK2L method) and the lower Runge-Kutta 
method of order one (the Euler-Cauchy or RK1L method)). The general formula for a lower Runge-Kutta 
method of order n is given next. 

Following Fehlberg (ref. 15), the integration formulas for a lower Runge-Kutta method of order n can 
be written as 


y(t + h) = y(t) + h^c K f K 

K=0 


(34) 


and 


f K =f 


K-l 

'+Vt,y + XlW* 

X=Q 


(35) 


The coefficients a K , PkX, and c K are determined by using Taylor’s expansion. The values of these coefficients 
for lower order Runge-Kutta methods of orders one and two, taken from Fehlberg (ref. 15), are listed in tables 
HI and IV. Complete derivation of these formulas and estimation of the coefficients may be found in this 
reference. 


Exponential Integration Method (or RK2EXP Method) 

Several explicit and implicit exponenticl integration algorithms have been proposed from time to time 
for the integration of the “stiff’ constitutive equations of viscoplastic models (refs. 11 to 13). These 
exponential integration algorithms are claimed to be considerably more stable (explicit algorithms) or 
unconditionally stable (implicit algorithms) in the “stiff’ region of integration provided that the integrand 
satisfies certain mathematical requirements. Applying some implicit exponential integration algorithms to 
stress analysis of structural engineering problems involving viscoplastic models (ref. 27) has shown them to be 
more efficient than the conventional integration methods. 

Explicit integration algorithms offer several advantages for the stress analysis of structural engineering 
problems involving many degrees of freedom. Therefore, an integration strategy based on an explicit 
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exponential integration method (ref. 16) was developed for application to structural (and life) analysis 
problems involving viscoplastic models. The exponential integration method developed in reference 16 for use 
in conjunction with the explicit Runge-Kutta method is shown to possess better stability characteristics than 
the conventional explicit integration methods. A built-in capability of the method automatically determines 
whether the problem is stiff or nonstiff and efficiently integrates it accordingly. This exponential integration 
algorithm has successfully been applied to several mathematically stiff problems and shown to be very 
efficient. A brief description and derivation of this integration method is followed by its numerical 
implementation and the description of the integration strategy for application to viscoplastic models. 

Rigorous mathematical details are avoided to keep the presentation simple. Readers interested in complete 
details may find these in reference 16. 

Integrating equation (32) from t to t+h yields 


t+h 

y(t + h) = y(t ) + J f[s, 

S=! 


(36) 


This is the fundamental theorem of integral calculus. 

To develop an exponential integration scheme, start with the Laplace integral 


A 

/ = J e~ xs F(s) ds , where x > 0 (37) 

0 


and F(s ) is of the exponential order. Write 


F(s) = F(s)e ms 


or 


F(s) = e~ 


' |f(0) + Jf'(O) + mF(0)]s + of* 2 )} 


by using Mclaurnri2s expansion. Choose m such that 

F(0) + mF(0) = 0 or m = -F'(0)/^(0) 


Therefore, 

A 

I = ^ e - (x ^ m)s F{s)e ms ds 

0 

Integrating gives 

x + m L J 

This result is used to accomplish the exponential integration of equation (32). 


(38) 


(39) 


(40) 


(41) 


(42) 
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Substituting s = t + u into equation (35) yields 

h 

y(t + h) = y(t)+ J e~ mu y' (t + u)e mu du (43) 

u=0 

Also 

e mu y'(t + u) = y'(t) + [y " (t) + my '(r)]« + O^w 2 j (44) 

Again, choose m such that 

y"(t)+my'(t) = 0 or m = -y"(t)/y'(t) (45) 

Using equations (42) to (45) gives 

1 - e~ mh 

y(t + h) = y(t)+y'(t) (46) 

m 

For m > 0, equation (46) yields an approximation to equation (32) that is useful for a large mora small h. For 
m < 0 and decaying exponential (stiff) behavior this method is not applicable. In this case use an explicit 
integration (e.g., second-order Runge-Kutta) method as the nonstiff integration method. It is important to note 
that by examining the sign of the exponent m the RK2EXP method automatically determines whether the 
given differential equation is stiff (m > 0) or nonstiff (m < 0) in the region of integration and integrates it 
accordingly. 

NUMERICAL IMPLEMENTATION 

The procedure for implementing the exponential integration method is given here in a step-by-step 

form: 

(1) Calculate derivatives at the start of the step. 

f[t,y(t)] = d M 

(2) Calculate the Euler values. 

z(t + h) = y(t) + h f[t,y(,t)] 

(3) Calculate the new derivatives. 

f[t + h,z(t + h)] = d new 

Note that these calculations are required also for the classical Runge-Kutta method. 

If <fold 4 s 0, 

Calculate m = {d 0 \& - d nev/ )l(h x d old ) 
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If m > 0, 

Calculate the exponential approximation 

y(t + h) = y(t) + d M x ( 1 - exp(-m x h)fm 


Endif 

Endif 


If ^old = 0 or m < 0, 

Calculate the RK2 approximation 

y{t+h)= MO + [hn] x (</ old + d DCV/ ) 


Endif 


ESTIMATION OF LOCAL ERROR (IMBEDDING TECHNIQUE) 

To estimate the local error in the solution, an “imbedding” technique is employed. In this technique 
the local error is defined as the difference in the solutions obtained by using two methods of different orders. 

For the RK2 method the local error is the difference between the RK2 and readily available Euler values. For 
the exponential method it is estimated by employing the exponential approximation (see ref. 16) 

y a (,t + h) = y(t) exp[y(fWy(0]+o(/» 1 2 j (47) 


where the subscript a refers to an approximation to the integration method. 


AUTOMATIC TIME-STEP INTEGRATION 

Two schemes were adopted to make the integration strategy self-adaptive, that is, capable of 
automatically picking suitable time steps that preserve the stability and accuracy of the solution: 


(1) Scheme I. Employed herein, scheme I defines the error in the solution during the current time 

step as 


_N 


* M 


(48) 


(2) Scheme D. It defines the error in the solution during the current time step as 


F 

H 


(49) 


Generalization to a system of n equations is easily obtained by defining an error norm e as 


e = 


n 

X e r, e r, 

i=i 


V /2 


(50) 
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The integration strategy runs as follows: 

(1) Prescribe the upper e u and lower limits for the error norm e. 

(2) If e > e u , halve the time step to hi 2. (If necessary, repeat this until a value of h is obtained that 
satisfies the condition e < e u .) 

(3) If e < ei, double the current time step to 2 h. 

(4) If ei < e < e u , retain the current time step. 

Proceed to the next step calculations with the time step rendered by the preceding and so on. 


APPLICATIONS 

To illustrate the applications of the self-adaptive time integration strategy in conjunction with the four 
explicit integration methods mentioned earlier, computer programs in FORTRAN were developed and 
applied to several uniaxial problems. The uniaxial problems considered consisted of tensile, relaxation, 
and cyclic thermaL/mechanical loadings. 

Figures 1(a), (b), and (c) exhibit tensile, relaxation, and cyclic loadings for the F-V model. The 
tensile loading (fig. 1(a)) consisted of a stress rate of 6.10X 10" 3 MPa/s at 149 °C. For relaxation loading 
(fig. 1(b)) the specimen was deformed to a strain of 0.015 in 1000 s and then allowed to relax at the same 
strain for 10 000 s. For in-phase cyclic thermomechanical loading (fig. 1(c)) the specimen was cycled 
between +0.00375 strain in 1000 s. The temperature of the specimen was cycled (in phase with the strain) 
between 200 °C (at -0.00375 strain) and 500 °C (at 0.00375 strain) in 1000 s. 

The tensile, relaxation, and cyclic thermal/mechanical loadings for the B-P model are shown in 
figures 2(a), (b), and (c), respectively. The tensile loading (fig. 2(a)) consisted of loading the specimen 
to 0.010 strain in 120 s at 871 °C. For relaxation loading (fig. 2(b)) the specimen was loaded to 0.008 
strain in approximately 100 s at 871 °C and then allowed to relax for 1000 s. For cyclic thermal/ 
mechanical loadings (fig. 2(c)) the strain and temperature were cycled in phase between ±0.004 and 538 
and 760 °C, respectively. The time for one cycle was 160 s. 


RESULTS AND DISCUSSION 

The constitutive equations of the F-V and B-P viscoplastic models were integrated in conjunction with 
the Runge-Kutta second-order method (RK2), the lower Runge-Kutta method of order one or the Euler- 
Cauchy (RK1L), the lower Runge-Kutta method of order two (RK2L), and the exponential (RK2EXP) 
integration method. The respective tensile, relaxation, and cyclic thermal/mechanical loadings were 
considered in numerical computations for the F-V and B-P models. The results of these computations are 
shown in figure 3 for the F-V model and in figure 4 for the B-P model. 

These computations were performed to ensure that the implementation of the self-adaptive integration 
strategy in conjunction with these four explicit integration methods yielded accurate results. The 
computations presented in figures 3 and 4 were carried out by assigning very small error tolerances for 
integration for different types of loadings (described earlier). The curves marked as “exact” in these 
figures were obtained by using the Runge-Kutta, second-order method (without self-adaptive integration 
strategy) and choosing a very small time step for integration of the constitutive differential equations of 
the two viscoplastic models. Figures 3 and 4 show that for both viscoplastic models and all uniaxial 
problems (tensile, relaxation, and cyclic loadings), the self-adaptive integration strategies produced 
results that are in close agreement with the “exact” results. This conclusion gives confidence in accurate 
implementation of the self-adaptive integration strategies for the four explicit integration methods. 

To test the efficiency of different numerical integration algorithms in conjunction with the self- 
adaptive time integration strategy described earlier, a large number of cases, involving the same loadings 
and viscoplastic models, were run. Different values and combinations of error tolerances were considered 
in numerical computations involving different types of loadings. Because it is not possible to list all the 
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results, the most significant results from these calculations for the three types of loadings are given in 
table V for the F-V model and in table VI for the B-P model, 

A criterion must be established for comparing the efficiencies and accuracies of different integration 
algorithms. The following criterion was adopted in the present study. First, an arbitrary and desired 
accuracy (or acceptable error) in the solution was assumed. Then, the error tolerances for each integration 
algorithm and error control scheme were varied until the required accuracy was achieved in the solution. 
This was done for each type of loading and viscoplastic model. The execution time (CPU time) for each 
case was recorded. The execution time for tracing a given loading path when using the RK2EXP (or any 
other) integration algorithm depended not only on the error tolerances but also on the scheme being used 
to define the time steps during the integration. Of the execution times thus obtained, the smallest time for 
each integration algorithm and the corresponding error tolerances, scheme, etc., were selected. These 
times are listed in tables V and VI for all loadings and viscoplastic models. These tables readily reveal the 
most efficient integration algorithm for a given viscoplastic model and given loading. 

It is advantageous to define the format of tables V and VI at this stage. The tables have identical 
formats. Column 1 lists the integration algorithms investigated. Column 2 indicates which error norm 
(scheme), equation (48) or equation (49), was used to control the time step for the self-adaptive time 
integration strategy. The prescribed upper error limits RTOL w for the Runge-Kutta methods and ETOL w for 
the exponential integration method are shown in columns 3 and 4, respectively. These values represent the 
value of e u (see section Automatic Time-Step Integration) for the two integration methods. An additional 
advantage of using the RK2EXP method, as may be noted from columns 3 and 4, is that it allows various 
combinations of error tolerances RTOL w and ETOL u , making the RK2EXP integration method more 
versatile than the other integration methods. Column 5 shows the total number of steps (iterations) 
required by a given method to trace the given loading path (cycle). The number of steps for the RK2EXP 
method in column 5 is split into two parts, a and b. The superscript a on a number denotes the number of 
times the RK2 method was used, and the superscript b denotes the number of times the exponential 
method was used. Column 6 lists the execution times (CPU times) in seconds (on a Cray-YMP computer) 
taken by various integration methods to trace a given loading path. Finally, the accuracy of a given 
integration method is indicated in column 7. This column lists the error in the result obtained by using a 
given integration algorithm. The error is defined as 


Error = (Y c - Y e )IY e 


where Y c denotes the solution obtained by using the current integration algorithm and Y e designates the 
“exact” solution obtained by using the RK2 method with a constant and small time step as described 
earlier. 


Results Using Freed- Verrilli Model 

Tensile loading . — The tensile loading shown in figure 1(a) and the F-V model were used to generate 
the results listed in table V(a). The RK1L and RK2L integration algorithms were found to be less efficient 
than the RK2 method. The execution times for the former two methods were 10.002 and 12.949 s, respec- 
tively, whereas for the same accuracy (0.460x10“^) the time taken by the RK2 method was 9.871 s. The 
most interesting result for the tensile loading case is obtained from the fourth row of the table. It shows the 
RK2EXP method to be significantly more efficient (by a factor of about 2) and slightly more accurate 

(0.416x 10”3) than the RK2 method (0.460x 10“3). Fewer iterations (40 352) were performed by the 
RK2EXP method to trace the tensile loading path than by RK2 and the other integration methods 
(104 091). 

Relaxation loading . — The results of integrating the F-V model over the relaxation loading path of 
figure 1(b) by the different integration methods are shown in table V(b). For the same accuracy the RK1L 
and RK2L methods were comparable to or less efficient than the RK2 method. The RK2EXP method, for a 
comparable accuracy, took less execution time (1.096 s) than the RK2 method (1.205 s). The RK2EXP 
method was about 10 percent more efficient (faster) than the RK2 method and took fewer steps (8381) 
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than the RK2 method (10 802). Using the RK2EXP method in this case is clearly advantageous as, for 
comparable accuracy, it is more efficient than the RK2 method. 

Cyclic loading .— Table V(c) summarizes results for the cycling loading (see fig. 1(c)) obtained with 
the F-V model. The RK2, RK1L, and RK2L integration methods gave almost the same accuracy (column 
7). However, RK1L and RK2L were not as efficient as the RK2 method. The execution times of 3.554 and 
4.487 s for the RK1L and RK2L methods were higher than that for the RK2 method (3.508 s). Comparing 
the results for the RK2 and RK2EXP methods (rows 1 and 5) shows that the RK2EXP method took less 
execution time (3.259 s) than the RK2 method (3.508 s). However, the error was larger (0.343x10-3) for 
RK2EXP than for RK2 (0.227x 10“3). Therefore, for cyclic loading the RK2 and RK2EXP methods are 
comparably advantageous and either may be used. 


Results Using Bodner-Partom Model 

Tensile loading . — The results for the tensile loading of figure 2(a) and the B-P model are shown in 
table VI(a). The total numbers of iterations performed using the RK2, RK1L, and RK2L integration 
methods were the same. However, the RK2EXP integration method took only 2129 steps to trace the given 
tensile loading path. The most important results of the tensile loading case are obtained from columns 6 
and 7. Column 6 shows that for the present case the RK2EXP method took significantly less execution 
time (0.223 s) than the RK2 (1.166 s), RK1L (1.183 s), and RK2L (1.620 s) integration methods. Also the 
RK2EXP method yielded a more accurate solution (error of 0.104xl0" 3 ) than the RK2 method (error of 
0.174xl0" 3 ). Therefore, using the RK2EXP method to integrate the constitutive equations of the B-P 
model is of significant advantage in this case. 

Relaxation loading . — Table VI(b) displays the results for relaxation loading (see fig. 2(b)) obtained by 
the different integration methods. For a given accuracy the execution times were higher for the RK1L and 
RK2L integration methods than for the RK2 method. The execution time was also marginally higher for 
the RK2EXP method (6.653 s) than for the RK2 method (6.270 s). The total numbers of iterations 
performed by the different integration methods to trace the same relaxation loading path were comparable. 
In this case the RK2EXP method was found (at best) to be comparable to the RK2 method. 

Cyclic loading . — The results of integrating the constitutive equations of the B-P model for different 
integration methods are listed in table VI(c) for the cyclic loading of figure 2(c). For a given accuracy 
the execution times were higher for the RK1L, RK2L, and RK2EXP integration methods than for the RK2 
integration method. The total number of iterations performed for the RK2, RK1L, and RK2L methods were 
the same (96 058). However, for the RK2EXP method, it was 96 059. The RK2EXP method does not offer 
any particular advantage with regard to efficiency and accuracy in this case. 


CONCLUSIONS 

The results presented for both viscoplastic models and three types of loadings show that, for 
comparable accuracy, the lower Runge-Kutta methods of orders one and two are less efficient than the 
Runge-Kutta, second-order and exponential integration methods. The lower Runge-Kutta methods are, 
therefore, not advantageous for application to viscoplastic models. The results also show that, for 
comparable accuracy, the efficiency of the integration algorithms depends on the type of loading. 

For example, for tensile loading, for comparable accuracy and for both viscoplastic models, the 
exponential integration method is much faster than the Runge-Kutta, second-order method. However, for 
the relaxation and cyclic loadings and for comparable accuracy, the exponential method is at best as 
efficient as the Runge-Kutta, second-order method. 

Because, in general, for both viscoplastic models and all the loadings considered, the exponential 
integration algorithm provides significantly higher (tensile loading) or comparable (relaxation and cyclic 
loadings) efficiency of any other integration algorithm, its use for integrating the “stiff’ equations of 
viscoplastic models may be advantageous. Note that the results presented in this paper pertain to uniaxial 
loadings and uniaxial forms of the viscoplastic models. Computer programs in FORTRAN were developed 
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to test the efficiency and accuracy of various integration algorithms with self-adaptive strategy. The 
structural engineering problems, however, involve multiaxial loadings and forms of the viscoplastic 
models. To judge the efficiency and accuracy of the exponential integration algorithm for application to 
these problems, it should be implemented into a finite element code. The process of implementing the 
exponential integration algorithm with self-adaptive strategy into a finite element code is in progress. The 
results of this work will be reported subsequently. 
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TABLE I.— FREED- VERRILLI MODEL 
CONSTANTS FOR COPPER 
[From reference 17.] 


Constant 

Unit 

Value 

E 

MPa 

165 000 - 125 T 

a 

°C-1 

15x10-6 + 5x10-9 T 


— 

0.34 


S " 1 

50x106 


MPa 

1.5 


MPa 

500 


MPa 

5000 

mm 

MPa 

25 exp(- 77300) 

n 

— 

5 

Q 

J/mole 

200 000 

s 

MPa 

14.3 

T m 

K 

1356 
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TABLE n.— BODNER-PARTOM MODEL CONSTANTS FOR B1900+Hf 

[From reference 18.] 

(a) Temperature-independent constants 


m\ =0.270 MPa -1 
m 2 = 1.520 MPa- 1 

ai =0.0 
Z { = 3000 MPa 
Z 3 = 1150 MPa 
r, = r 2 = 2 
D 0 =lxl0 4 s-' 


(b) Temperature-dependent constants 


Constant 

Temperature, °C 

<760 

871 

982 

1093 

n 

Zo (MPa) 

A[ =A 2 (S-*) 

Z 2 = Z 0 (MPa) 

1.055 

2700 

0 

2700 

1.03 

2400 

0.0055 

2400 


0.7000 

1200 

0.2500 

1200 


(c) Elastic constants (MPa; T in °C) 


E= 1.987x10 s + 16.78 T - 0.1034 t 2 + 1.143xl0“ s T 3 
G = 8.650x1 0 4 - 17.58 T + 0.2321 T 2 - 3.464X 10" 5 T 3 


TABLE HI. — COEFFICIENTS FOR LOWER RUNGE-KUTTA 
METHOD OF ORDER ONE (EULER-CAUCHY 
OR RK1L METHOD) 


K 

“k 

PkX 

(X=0) 

m 

0 

0.0 


0.0625000 


i 

.5333333 

0.5333333 

.5333333 


TABLE IV — COEFFICIENTS FOR LOWER RUNGE-KUTTA 
METHOD OF ORDER TWO (RK2L METHOD) 


K 

<*k 

PkX 


X 

0 

1 

0 

0.0 

0.0 


0.2401795 

i 

.2500000 

.2500000 


.0303030 

2 

.6750000 

-.2362500 

0.9112500 

.7295174 
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TABLE V.— RESULTS FOR FREED- VERRILLI MODEL 


(a) Tensile loading 


Method 

Scheme 

Prescribed 
upper error in 
Runge-Kutta 
method, 
RTOL, 

Prescribed 
upper error in 
exponential 
method, 
ETOL u 

Iterations 

Execution 
(CPU) time, 
s 

Error 

RK2 

i 

10-3 


104 091 

9.871 

0.460x10-3 

RK1L 

i 

10-3 


104 091 

10.002 

.460x10-3 

RK2L 

i 

10-3 


104 091 

12.949 

.460x10-3 

RK2EXP 

n 

10-7 

5x10-6 

40 352 (34743+36 878b) 

4.638 

.416x10-3 


(b) Relaxation loading 


RK2 

I 

mm 


10 802 

mm 

0.352x10-3 

RK1L 

I 

MM 


10 802 

mam 

.352x10-3 

RK2L 

I 

H 


10 802 


.352x10-3 

RK2EXP 

n 

MM 

io - 4 

8381 (10013+7380^) 

1.096 

.364x10-3 


(c) Cyclic loading 


RK2 

I 

10-2 


32 105 

3.508 

0.227x10-3 

RK1L 

I 

10-2 


32 102 

3.554 

.223x10-3 

RK2L 

I 

10-2 


32 105 

4.487 

.223x10-3 

RK2EXP 

n 

10-6 

io- 4 

26 160 (11 1573+15 003b) 

3.259 

.343x10-3 


a Number of times Runge-Kutta method was used. 
^Number of times exponential method was used. 
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TABLE VI. — RESULTS FOR BODNER-PARTOM MODEL 


(a) Tensile loading 


Method 

Scheme 

Prescribed 
upper error in 
Runge-Kutta 
method, 
RTOL// 

Prescribed 
upper error in 
exponential 
method, 
ETOL w 

Iterations 

Execution 
(CPU) time, 
s 

Error 

RK2 

i 

10-2 


11 510 

1.166 

0.174x10-3 

RK1L 

i 

10-2 


11 510 

1.183 

.174x10-3 

RK2L 

i 

10-2 


11 510 

1.620 

.174x10-3 

RK2EXP 

n 

io-i 

10-1 

2129 (1948 a +181 b ) 

.223 

.104x10-3 


(b) Relaxation loading 


RK2 

I 

10-3 


61 000 

6.270 

0.525x1 0 -4 

RK1L 

I 

10-3 


61 000 

6.299 

.525X10- 4 

RK2L 

I 

10-3 


61 000 

8.712 

.525x10-4 

RK2EXP 

D 

10-3 

10-3 

61 008 (50 107a+10 901b) 

6.653 

.525x10-4 


(c) Cyclic loading 


RK2 

I 

10-5 


96 058 

■SB 

0.858x10-5 

RK1L 

I 

10-5 


96 058 

1 

.858x10-5 

RK2L 

I 

10-5 


96 058 

mmm 

.858x10-5 

RK2EXP 

n 

10-5 

10“ 5 

96 059 (62 941®+33 118b) 

■SB 

.858x10-5 


a Number of times Runge-Kutta method was used. 
^Number of times exponential method was used. 
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Strain 


Figure 3. — Stress-strain and relaxation curves and 
hysteresis loops for Freed-Verrilii model (copper). 

(a) Stress-strain curve. Stress rate, 6.1 0x1 0 -3 MPa/s; 
temperature, 149 °C. (b) Relaxation curve. Strain 
rate, 1.50x1 0-^/s; temperature, 316 °C. (c) Hysteresis 
loops. Strain rate, 1.50x1 O^/s; strain range, ±0.00375; 
temperature, 200 <=> 500 °C. 





Figure 4. — Stress-strain and relaxation curves and hysteresis 
loops for Bodner-Partom model (B1 900+ Hf). (a) Stress-strain 
curve. Strain rate, 8.333x1 0*^/s; temperature, 871 °C. 

(b) Relaxation curve. Strain rate, 8.333x1 O^/s; maximum 
strain, 0.008; temperature, 871 °C. (c) Hysteresis loops. Strain 
rate, 1.0x1 0"4/s; strain range, ±0.004; temperature, 538 « 
760 °C. 
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